from __future__ import print_function
import SimpleITK as sitk
import numpy as np
import os
import matplotlib.pyplot as plt
#matplotlib inline

output_path = "imageOutputPath"

#########################################################################################
#########################################################################################
### 查看slice

itk_img = sitk.ReadImage('mhdDir')

img_array = sitk.GetArrayFromImage(itk_img)   # indexes are z,y,x (notice the ordering)
num_z, height, width = img_array.shape        # heightXwidth constitute the transverse plane

#MHD值的坐标体系是体素，以mm为单位（dicom的值是GV灰度值）
#结节的位置是CT scanner坐标轴里面相对原点的mm值，需要将其转换到真实坐标轴位置
#可以使用SimpleITK包中的 GetOrigin() GetSpacing()
origin = np.array(itk_img.GetOrigin())      # x,y,z  Origin in world coordinates (mm)
spacing = np.array(itk_img.GetSpacing())    # spacing of voxels in world coor. (mm)

# 以z轴=10切片显示
plt.imshow(sitk.GetArrayViewFromImage(itk_img)[10,:,:], cmap=plt.cm.Greys_r)
plt.axis('off')

#########################################################################################
#########################################################################################

# 分割出的结果切片存成bmp
# 原始方法跟3*3*3卷积核方法结果比较

from PIL import Image

itk_img_o = sitk.ReadImage('originalMehtodResult.mhd')
itk_img_3 = sitk.ReadImage('3^3MethodResutl.mhd')

img_array_o = sitk.GetArrayFromImage(itk_img_o) 
img_array_3 = sitk.GetArrayFromImage(itk_img_3) 

num_o, height_o, width_o = img_array_o.shape
num_3, height_3, width_3 = img_array_3.shape
img_array_m = np.zeros([num_o,height_o,width_o])

for zz in range(0, num_o-1):
        for hh in range(0, height_o-1):
            for ww in range(0, width_o-1):
                img_array_m[zz,hh,ww]=(img_array_o[zz,hh,ww]*128 + img_array_3[zz,hh,ww]*64)

def Matrix2Image(data):
    #data = data*255
    new_img = Image.fromarray(data.astype(np.uint8))
    return new_img

def Matrix2Imageo(data):
    data = data*128
    new_img = Image.fromarray(data.astype(np.uint8))
    return new_img

def Matrix2Image3(data):
    data = data*64
    new_img = Image.fromarray(data.astype(np.uint8))
    return new_img

data = img_array_m[zz/2,:,:]
Matrix2Imageo(img_array_o[zz/2,:,:]).save('o.bmp')
Matrix2Image3(img_array_3[zz/2,:,:]).save('3.bmp')
Matrix2Image(data).save('cool.bmp')

#########################################################################################
#########################################################################################

# 切片 输出npy
# 切片可输出array为npy文件
# SimpleITK的img_array的数组不是直接的像素值，而是相对于CT扫描中原点位置的差值，需要做进一步转换
# python 中图像用SimpleITK和numpy.ndarray表示的差异
# http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/Python_html/03_Image_Details.html
# 主要是坐标轴顺序 直接换一下就好

b = sitk.GetArrayViewFromImage(itk_img)[23,:,:]
np.save(os.path.join(output_path,"images_23.npy"),b)

#########################################################################################
#########################################################################################

# 评分
# 逐点计算

itk_img_0 = sitk.ReadImage('originalMehtodResult.mhd')
itk_img_1 = sitk.ReadImage('3^3MethodResutl.mhd')

img_array_0 = sitk.GetArrayFromImage(itk_img_0) 
img_array_1 = sitk.GetArrayFromImage(itk_img_1) 

num_z_0, height_0, width_0 = img_array_0.shape       
num_z_1, height_1, width_1 = img_array_1.shape

if (num_z_0 != num_z_1)or(height_0 != height_1)or(width_0 != width_1):
    raise ValueError
else:
    v1=0 
    v2=0
    v3=0
    for zz in range(0,numz_0-1):
        for hh in range(0,height_0-1):
            for ww in range(0,width_0-1):
                v1= v1+ img_array_0[zz,hh,ww]
                v2= v2+ img_array_1[zz,hh,ww]
                v3= v3+ img_array_0[zz,hh,ww]* img_array_1[zz,hh,ww]
    print('gt:', v1)
    print('test:', v2)
    print('overlap:', v3)
    print('dice:',v3*2.0/ (v1+ v2))
a = v3*2.0/ (v1+ v2)
print a 





